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Abstract 

Lattice models are central to the physics of ultracold atoms and condensed matter. Generally, 
lattice models contain time-independent hopping and interaction parameters that are derived from 
the Wannier functions of the noninteracting problem. Here, we present a new concept based 
on time- dependent Wannier functions and the variational principle that leads to optimal time- 
dependent lattice models. As an application, we use the Bose-Hubbard model with time-dependent 
Wannier functions to study an interaction quench scenario involving higher bands. We find a 
separation of time-scales in the dynamics. The results are compared with numerically exact results 
of the time-dependent many-body Schrodinger equation. We thereby show that - under some 
circumstances - the multi-band nonequilibrium dynamics of a quantum system can be obtained 
essentially at the cost of a single-band model. 
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Quantum phenomena of particles in periodic and quasi-periodic potentials are central 
themes in theoretical physics. Alone the question about the nature of the many-body ground 
state in such potentials has been investigated for decades [lj]. Ultracold atoms in optical 
lattices are highly controllable realizations of such many-body systems and allow to directly 
monitor their nonequilibrium dynamics. Moreover, they promise to provide insight into the 
physics of real solids. Lattice models play a crucial role in the current physical understanding 
of such systems. At the heart of a lattice model is the idea of lattice site localized orbitals 
which are commonly referred to as Wannier functions 

a a. 

In a lattice model Wannier 



amiltonian. Among the 
4J and fermionic [5!] ver- 



functions are used as a single-particle basis for the many-body H; 
most famous lattice models are certainly the single-band bosonic 
sions of the Hubbard model, both of which predict phase transitions between insulating and 
superffuid phases. Recently there has been large interest in the multi-band physics of ul- 
tracold atoms by both theorists and experimentalists a 
already couple the ground state to higher bands 



ike 



15) . Interparticle interactions 



and it was shown recent 



y that 



even very weak interactions lead to multi-band physics in nonequilibrium problems 141 ] . 

Multi-band phenomena seem to be abundant and we would like to address the question 
of the optimal way to deal with them theoretically. This calls for optimal lattice models 
that incorporate multi-band physics effectively. Here, we present a completely new concept 
that leads to variationally optimal lattice models for multi-band nonequilibrium dynamics 
and statics problems. Our idea is to allow the Wannier functions of a given lattice model 
to depend on time and to determine the many-body dynamics from the variational prin- 
ciple. The concept is applicable to bosons and fermions as well as any number of lattice 
sites, particles and bands. As an application, we use the Bose-Hubbard model with time- 
dependent Wannier functions to study a sudden quench of the interaction in a double-well 
problem that involves higher bands. We find a separation of time-scales in the dynamics. 
On short time-scales higher-band excitations lead to rapid density oscillations, not captured 
on the single-band level. On longer time-scales we observe oscillations between condensed 
and fragmented states. The results are compared with numerically exact results of the 
time-dependent many-body Schrodinger equation. We thereby show that when the lattice 
is sufficiently deep and the interaction couples to higher bands the multi-band nonequilib- 
rium dynamics of a quantum system can be obtained essentially at the cost of a single-band 
model. Already the physics of this example proves to be much richer than what the single- 
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band Bose-Hubbard (BH) model can access. 

Consider a one-dimensional (ID) lattice of M sites. We make the following ansatz ior the 
time- dependent Wannier functions 

w k (x, t) = Y, d a k(t)wt(x), k = 1, . . . , M, (1) 

where w k (x) denotes the conventional Wannier function of the band a at lattice site k and 
the d k (t) its time-dependent amplitude. If v is restricted to one, the ansatz ([[]) reduces to 
the conventional Wannier functions of the lowest band. Using (JTj) as a basis, the ansatz for 
the many-boson wave function becomes = Yln^n{t) \n',t), where the sum is over all 

permanents \n;t) = \ni, ri2, ■ ■ ■ , rijwj t) of N bosons residing in M time-dependent Wannier 
functions. We write b k (t) (b k (t)) for the operator that creates (annihilates) a boson in 
the time-dependent Wannier function Wk(x,t) and n k (t) = b\.(t)b k (t) for the corresponding 
occupation number operator. The bosonic commutation relations [bk{t),tf {t)] = 5k q hold at 
all times. For simplicity we use dimensionless units in which h = m = 1. Analogous to the 
derivation of the BH model, the Hamiltonian of the time-dependent BH (TDBH) model 

Htdbh = J2\ - Jkk+i(t)P k (t)h+i(t) - J* kh+ i(t)bl+i(-t)h(t) 

k=l ^ 

+e k {t)h k {t) + [h k (t) {n k {t) - I)]) (2) 

Xi X j ) . 



can then be derived from the full many-body Hamiltonian YuiLi h( x i) + -^o J2f<j 
Here, h(x) = — + V(x) is the one-body part of the Hamiltonian including the 
external potential V(x). The time-dependent parameters in Eq. §2§ are defined as 
Jkk+iit) = - J dxw k (x,t)h(x)w k+ i(x,t), e k (t) = j dxwl(x,t)h(x)w k (x,t) and U kkkk {t) = 
Ao / dx\w k (x, t)\ 4 . Contrary to the BH model, the hopping parameter J kk+ i(t), the one- 
body energy e k (t) and the interaction parameter U kkkk (t) are now time-dependent and can 
vary from one lattice site to another. Note also that the hopping parameter J kk+ i(t) will 
generally not be real. It is then possible to derive coupled equations of motion for the time- 
dependent Wannier functions and the coefficients of the many-body wave function from the 



variational principle 16l. Il7fl: 



i\w k (t)) = P k {t) 



Mw k (t)}+ V ^h\ Wl (t)} + ^^U kk (t)\w k (t)} 
^ Pkk{t) pkk{t) 



l=k±l 



iC(t) = H(t)C(t), (3) 



where C(t) = {Cait)}. The derivation and a discussion of Eqs. (j3J) are given in the ap- 
pendix. The operators P k (t) = Y^b=i \ w k)( w k I ~~ \ w k(t)) (wk{t)\ appearing in Eq. (j3J) are 
projection operators, Ukk(x,t) = Xo\wk(x, t)\ 2 and H^>(£) = (n\ t Htdbh n';t^. Similarly, 
the quantities p kl (t) = (*(t)|SJ(t)S,(t)|*(t)) and p kkkk (t) = (*(t)|SJ(t)St(t)S fc (t)S fc (t)|*(t)> 
are matrix elements of the first- and second-order reduced density matrices. We stress that 
the optimal lattice model governed by Eqs. ([U13]) preserves all conservation laws. In par- 
ticular, particle-number is a good quantum number, energy is conserved when the original 
Hamiltonian is time- independent, and spatial symmetries of this Hamiltonian are preserved 
in time for initial conditions preserving them. Summarizing the above, we have obtained 
the TDBH model for which the wave function and the parameters of the Hamiltonian (j2J) 
evolve according to the variational principle. The dynamics is determined by the coupled 
equations of motion (J3J) . 

Now we would like to demonstrate the above ideas in an apparently simple two-site 
example. The dynamics of bosons on two lattice sites has recently been studied intensively 
using the BH model 18M28I]. A question that arises naturally in dynamical problems is the 
response of a quantum system to a perturbation. When the system is in the ground state and 
the perturbation is realized as a variation of one of the parameters in the Hamiltonian, the 
problem is known as a quantum quench. Quenches have recently received growing attention 



in the context of quantum gases in lattices 



29l432j. Generally, the quench dynamics is 



studied using single-band Hubbard models or simplifications thereof. Here, we study a 
quench going beyond the BH model and would like to investigate what new physics appears. 
Here, V(x) = Vocos 2 (a;) is an external potential on an interval of length 2ir that realizes 
a ring lattice of M = 2 sites (denoted L and R) with Vq — 25E r , where E r = 1/2 is the 
recoil energy. The splitting 2J = 2.08 x 10~ 3 of the ground and the first excited single- 
particle state of V(x) determines the Rabi oscillation period = tt/J. To characterize 
the interaction strength we use the parameter A = Xq(N — 1). The quench is then realized 
as a sudden change from A = to A = 0.6 in a system of N = 20 bosons, initially prepared 
in the noninteracting ground state of V(x). Within the BH model this quench corresponds 
to a change from U/ J — to U/ J — 25.8. We will now study and compare the dynamics of 
the TDBH and the BH model. For the TDBH model the number of bands v used for the 
expansion of each time-dependent Wannier function in Eq. ([1]) is increased until convergence 
is reached. In the present example we found that the results for J(t), U(t), the one-particle 
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density and the momentum distribution do not change visibly anymore for v ~ 5, and 
v — 10 was used throughout for full convergence. Furthermore, the spatial symmetry of the 
problem implies that only bands with even Wannier functions, a = 1, 3, 5, 7, 9, contribute. 
Moreover, J{t) remains real, eiit) = €n(t) = e(t), and Ullll{£) = U^^ii) = U(t) at all 
times. All TDBH computations were carried out on a standard desktop computer and none 
of the computations took longer than a few hours. 

On physical grounds it can be anticipated that the sudden increase of the (repulsive) 
interaction strength will lead to breathing oscillations of the one-particle density p^(x\x; t) 
in each of the wells. However, as the final interaction strength is not very strong these 
density oscillations should have a small amplitude only. It is important to note that for the 
BH model the ratio of the interaction and the hopping parameter U j J is the only relevant 
parameter and also that it is constant in time after the quench. In the TDBH model on the 
other hand we expect U(t)/J(t) to vary, because the Wannier functions can change in time. 
This is indeed the case. Fig. [T] shows the parameters U(t) and J(t) of the TDBH model and 
their ratio U(t)/J(t) as a function of time. As expected, right after the quench the density 
in each well expands and hence U(t) decreases, while J(t) increases. Both parameters then 
oscillate at a high frequency. Note the short time-scale. The interaction broadens the density 
and therefore U(t) {J{t)) is always smaller (larger) than its value at time t = 0. As J(t) 
depends on the overlap of adjacent time-dependent Wannier functions, it is very sensitive 
to any variation of their tails. It varies over almost 25%! The interaction parameter U(t) is 
sensitive to variations of the bulk density and varies only over about 4%. As a result, the 
ratio U(t)/J(t) of the TDBH model varies between 25.8 and 20, a range of about 25%. In 
contrast, for the BH model the ratio U/J = 25.8 is constant at all times as mentioned above. 
The rapid oscillations of U(t) and J(t) are only possible because the Wannier functions can 
evolve in the space of multiple bands and thus represent multi-band excitations. This clearly 
demonstrates the advantages of the TDBH model over the BH model in reproducing the 
physics of the quench studied here. 

We now turn to the question about the nature of the quantum state after the quench and 
to longer time-scales. We first focus on the eigenvalues of the first-order reduced density 
matrix, the natural occupation numbers (i), which determine the degree to which the 



system is condensed or fragmented 



33 



34J. In Fig. [2] (top) the natural occupation numbers 



of the TDBH model are shown together with those of the BH model. For comparison also 
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the numerically exact dynamics of the time-dependent many-body Schrodinger equation is 
shown, obtained using the multiconfigurational time-dependent Hartree for bosons method. 
For details see the literature 



17 



351 ] . Note that the time-scale is now more than a hundred 



times larger compared to Fig. [TJ It is clearly seen that all three results display oscilla- 
tions between fragmented and partially condensed states. For the first few oscillations the 
three results essentially coincide, which implies that for this problem the natural occupation 
numbers are not very sensitive to the previously discussed multi-band excitations. However, 
after a few oscillations between condensed and fragmented states the BH result deviates 
substantially from the exact one. On the other hand, the TDBH result follows the exact one 
closely for many oscillations. In Fig. [2] (bottom) the accumulated error of the BH model, 
defined as the accumulated difference — J* dt'\n^ BH {t') — n^ exact (t')\ between the largest 
natural occupation number of the BH model and the numerically exact result are shown, 
together with the analogous quantity of the TDBH model. The accumulated error of the 
TDBH model is always smaller than that of the BH model and grows more slowly, which 
means that also on longer time-scales the TDBH model captures the true physics much 
better than the BH model. 

We will now discuss the implications of the differences of the natural occupation numbers 
on observables. The rapid density oscillations discussed earlier have a very small amplitude, 
and consequently hardly any dynamics is visible in real space. However, a substantial 
dynamics occurs in momentum space. Fig. [3] shows the one-particle momentum distribution 
pW(k\k; t) of the exact, the TDBH and the BH result at t = where they all coincide, and 
at later times where differences have developed. The momentum distribution of the TDBH 
model is always closer to the exact one and follows it for a long time, whereas that of the 
BH model deviates quickly. Here, the TDBH model allows to obtain the precise multi-band 
nonequilibrium dynamics essentially at the cost of a single-band model. Of course, in order 
to obtain the exact result the terms that were neglected so far in Eq. (j2j) have to be included 
into the model, e.g., the terms responsible for correlated tunneling. Thus, the advantages 
of the TDBH model over the BH model have been clearly proven. 

In this work we have generalized the concept of lattice models by letting their Wannier 
functions become time-dependent. For such lattice models equations of motion can be de- 
rived from the variational principle, leading to variationally optimal lattice models that can 
efficiently incorporate multi-band physics. The concept is applicable to both nonequilibrium 
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dynamics as well as statics. In the latter case the ground state can be computed using imagi- 
nary time-propagation. The variational principle ensures that for identical initial conditions 
the optimal lattice model improves on the original lattice model, if the latter contains the 
terms of the full many-body Hamiltonian relevant to the problem at hand. Since the nu- 
merical effort of the optimal and standard lattice models are comparable, the former can 
also be used to assess the quality of the latter: if the respective results differ noticeably, the 
standard lattice model is inapplicable. As an application, we have presented the equations 
of motion for the TDBH model, i.e., the Bose- Hubbard model with time-dependent Wannier 
functions. We then studied a quantum quench problem in a double-well using the TDBH 
model. Multi-band excitations were found that resulted in high-frequency, large amplitude 
oscillations of the hopping and interaction parameters. Such phenomena are not accessible 
to the BH model. On a longer time-scale the multi-band excitations also affect the nature 
of the quantum state and the momentum distribution. We have compared the results of the 
TDBH and the BH model for the natural occupation numbers and the momentum distribu- 
tions with numerically exact ones of the many-body Schrodinger equation. For this problem, 
we find that the TDBH results follow the exact ones for a long time, while the BH results 
deviate quickly. Interestingly, the numerical effort for solving the equations of motion of the 
TDBH model is essentially the same as for the BH model, see the appendix for details. 

Summarizing, standard lattice models can be greatly extended at little extra cost through 
the use of time-dependent Wannier functions. It will be interesting to examine what other 
physical multi-band phenomena are accessible to such time-dependent lattice models, e.g., 
in disordered media or in the case of the fermionic Hubbard model. As a further direction 
we suggest to combine ID optimal lattice models with the recently developed adaptive time- 



dependent density-matrix renormalization group methods 



36M38I] to treat non-equilibrium 



dynamics of larger lattice systems, as is presently done utilizing standard lattice models. 

Appendix: Derivation and discussion of the equations of motion of the TDBH model 

We will now derive Eqs. ([3]) and discuss their properties. Using the ansatz for the many- 
body wave function = XlnCn(^) ^) given above, the action functional of the TDBH 
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model reads 

S[{Cn},{d a k }} 



dt< < * 



Htdbh — i 



d_ 
dt 



^-X)A«i(*)[<tf*WMt)>-l]}, (A.l) 



where the Lagrange multipliers are introduced to ensure orthonormality between the time- 
dependent Wannier functions. Note that Eq. (CQ) implies that time-dependent Wannier 
functions on different lattice sites are orthogonal by construction. In order to take functional 
derivatives it is useful to express the expectation value {^\H T dbh — *^|^) as an explicit 
function of the amplitudes {d k (t)} and expansion coefficients {Ca(t)}. On the one hand 
(^\Htdbh ~ i§i\^) can ^ e written as: 



Htdbh — i 



d_ 

dt 



Htdbh — i 



d_ 

dt 



.dC H {t) 



dt 



• (A.2) 



Using Eq. (1A.2|) we now require stationarity of the action functional with respect to varia- 
tions of the coefficients, = s ^»^ , and obtain: 



n: t 



.0 

Htdbh — iw; 

dt 



n';t)C r7 ,(t)=z 



dC H (t) 
dt ' 



(A.3) 



On the other hand {^>\Htdbh — i§f\^) can a l so be written as: 



ft i 9 

n-TDBH — 1-^7 

at 



*)= E EM^d»*(*)(^-<Mo^4)df(*) 

' k,q=k,k±l a,P ^ ' 

1 dC*(t) 

+9 E E a**® d T(t)4(t) u kak , k , k s di(t)d s k (t) - 1 J2 cdt)—g^, 



(A.4) 



k a/3-yS 



where h kaq /3 = (w k \h\w^), U kak p kJk 5 — X j dxw k {x)*w k {x)*w 5 k {x)w 5 k {x) and we have used 
(w k \wP) = SkqSa/3. Note that the matrix elements of the first- and second-order reduced 
density matrices are functions of the coefficients {Ca(t)} only. Using Eq. (1A.4jl we now 
require stationarity of the action functional with respect to variations of the time- dependent 
amplitudes, = Sd i^^ , and obtain 

E /°%(*)( w fcH w ?(0) +Pkkkk(t)(w k \U kk {t)\w k (t)) = p k (t)(w k \w k (t))+ip kk (t)(w k \w k (t)), 

q=k,k±l 

(A.5) 

for a — 1, . . . , v and k — 1, . . . , M, where we have used d k (t) = (w k \w k (t)) . The Lagrange 
multipliers fi k (t) can be determined from Eqs. (1A.5I) . which upon resubstituting the result, 
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multiplying each equation by \w%) and summing over a gives 



ipkk(t)P k (t)\w k (t)) = h(t) 



22 Pkq(t)h\w q (t)) + Pkkkk{t)Ukk{t)\w k (t)) 
_q=k,k±l 



(A.6) 



By including the additional phase factor e~ $ dt '( w k(t)\w k (t)) ^ Q ^ e definition of \wk{t)), we 
find that (wk{t)\wk(t)) = 0, and that Eqs. (IA.6[) and flA.3[) reduce to Eqs. 

Let us briefly discuss the numerical advantage of using time-dependent Wannier functions 
over conventional ones. The fundamental difficulty in solving quantum lattice models lies 
in the fact that the dimension of the many-particle Hilbert space grows very quickly as 
a function of the number of particles N and the number of lattice sites M. In the case 
of the BH model the length of the vector C is ( N+ n~ 1 )- This is also the length of the 
vector C in the TDBH model since we have only replaced conventional Wannier functions 
by time-dependent ones. However, since each of the time-dependent Wannier functions 
can take on any shape allowed by Eq. ([I]), a much larger Hilbert space is available to 
the variational principle. Note that the number of bands v in which each time-dependent 
Wannier function is expanded does not enter the above dimension formula. If v is restricted 
to one, the time-dependent Wannier functions become time-independent and only the second 
of Eqs. ([3]) needs to be solved. In this case Eqs. (j3J) reduce to the equations of motion of 
the BH model. If v > 1 the coupled system, Eqs. (E]), needs to be solved. Then, the 
additional numerical cost consists of propagating not only the coefficients {C^(t)}, but also 
the time-dependent Wannier functions {wk{x,t)}, i.e., M additional nonlinear equations 
have to be solved. For all but the smallest numbers of particles and lattice sites, this 
extra effort is much smaller than the effort needed to propagate the second of Eqs. ([3]). 
Note that the alternative of enlarging the Hilbert space by using multiple, say k > 1, 
bands of conventional (static) Wannier functions, would increase the length of the vector 
C to ( N+K ^ 1 ) • This illustrates why it is crucial to use time-dependent Wannier functions: 
time- dependent Wannier functions allow to keep the numerical effort of a time- dependent 
lattice model very close to that of the respective standard lattice model while efficiently 
incorporating multi-band physics. Since all parameters, i.e., the coefficients {Cg(i)} and the 
Wannier functions {wk(x,t)} are at all times determined by the variational principle, the 
lattice model is optimal. 
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Figure 1: (color online). Bose-Hubbard model with time-dependent Wannier functions (TDBH 
model). Shown is the time-dependence of the hopping and the interaction parameters J(t) and 
U{t) for a quench. Initially a system of N = 20 bosons is prepared in the noninteracting ground 
state of a two-site ring lattice, see text for details. At t = the interaction strength is switched on 
to £/(0)/J(0) = 25.8 and the resulting dynamics is monitored. Top: U(t)/J(t) oscillates rapidly 
with a large amplitude between about 25.8 and 20. Note the time-scale. U (t) / J(t) is always smaller 
than its initial value which is the value of the BH model (dashed line). Bottom: J(t) and U(t) 
oscillate in time. J it) is always larger, U(t) always smaller than its initial value. The oscillations 
are due to multi-band excitations induced by the quench. J(t) varies over almost 25%, U(t) over 
about 4% of its initial value. The hopping process is obviously much more involved than what 
the time-independent parameter J of the BH model can capture. The respective BH results are 
constant in time and equal to the values at t = 0. All quantities shown are dimensionless. 
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Figure 2: (color online). Quench dynamics in a double well as in Fig. [TJ Top: shown are the 
natural occupations of the BH (solid magenta lines) model, the TDBH (solid black lines) model 
and the exact (dashed blue lines) result as a function of time. Starting from a fully condensed state 
the dynamics is intricate and shows oscillations between partially condensed and fragmented states. 
The TDBH result follows the exact one closely and in particular reproduces the frequencies of the 
oscillations well. The BH result deviates quickly from the exact one. Bottom: the accumulated 
error (see text) of the TDBH model is always smaller than that of the BH model and grows much 
more slowly. All quantities shown are dimensionless. 
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Figure 3: (color online). Quench dynamics in a double well as in Fig. [TJ Shown is the one-particle 
momentum distribution of the BH (solid magenta lines) model, the TDBH (solid black lines) and 
the exact result (solid blue lines) at different times. The initial state is fully coherent and for a 
short time the three results have similar momentum distributions (not shown). The TDBH result 
is always much closer to the exact one than the BH result. Differences of the TDBH (BH) result 
relative to the exact one always occur when the natural occupations differ, see Fig. [2J All quantities 
shown are dimensionless. 
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